#include "yocto/math/fcn/zfind.hpp"
#include "yocto/ios/ocstream.hpp"

using namespace yocto;
using namespace math;

static const double K_AE2=10e-3;
static const double K_Henry = 29.41; //  -- atm/(mol/L)
static const double K1      = 4.45e-7;  
static const double K2      = 5.6e-11;
static const double K1p     = K1/K_Henry;

static const double L0 = 1000;
static const double Kr = 1.8e-8;
static const double Kt = 3.6e-6;
static const double C  = Kr/Kt;

typedef numeric<double>::function Function;

const double P0 = 40;

class Params
{
public:
    Params( double pH0 ) :
    P( P0 ),
    sig_AE2( this, & Params::get_sig_AE2),
    sig_NHE( this, & Params::get_sig_NHE),
    sig_Diff(this, &Params::regulate),
    sig_DelP(this, &Params::regulate2),
    h0( pow(10,-pH0) ),
    solver(),
    mu(0)
    {
        solver.xerr = 1e-5;
        Function F(this,&Params::get_diff);
        double r=1;
        while( F(r) <=0 ) r += 0.1;
        (double&)mu = solver(F,0,r);
    }
    
    double        P;
    Function      sig_AE2;
    Function      sig_NHE;
    Function      sig_Diff;
    Function      sig_DelP;
    double        h0;
    zfind<double> solver;
    const double  mu;
    
    
    
private:
    YOCTO_DISABLE_COPY_AND_ASSIGN(Params);
    
    inline double get_sig_AE2( double h )
    {
        const double b = K1p * (P/760.0) / h; // bicarbonates
        return b / (K_AE2+b);
    }
    
    inline double get_sig_NHE( double h )
    {
        const double x    = h/Kr;
        const double xp1  = x+1;
        const double Cx   = C*x;
        const double Cxp1 = Cx+1;
        const double num  = x*xp1 + L0*Cx*Cxp1;
        const double den  = L0 * (Cxp1*Cxp1) + (xp1*xp1);
        return num/den;
    }
    
    inline double get_diff( double r  )
    {
        return r * sig_AE2(h0) - sig_NHE(h0);
    }
    
    inline double regulate( double h )
    {
        return mu * sig_AE2(h) - sig_NHE(h);
    }
    
    inline double regulate2( double PiCO2 )
    {
        P = PiCO2;
        return mu * sig_AE2(h0) - sig_NHE(h0);
    }
    
};

int main()
{
    
    try
    {
        //----------------------------------------------------------------------
        // initialize AE2/NHE ratio 
        //----------------------------------------------------------------------
        double pH0  = 7.2;
        Params param(pH0);
        
        std::cerr << "-- mu= " << param.mu << std::endl;
        
        {
            ios::ocstream fp( "sig_diff.dat", false );
            for( double pH = pH0-1.5; pH <= pH0+1.5; pH += 0.01 )
            {
                const double h = pow(10,-pH);
                fp("%g %g\n", pH, param.sig_Diff(h) );
            }
            
        }
        
        //----------------------------------------------------------------------
        // set pH to extrema pH
        //----------------------------------------------------------------------
        double pHmin = 7;
        param.h0 = pow(10,-pHmin);
        
        {
            ios::ocstream fp( "sig_delp.dat", false );
            for( double P = P0; P <= 2*P0; P += 0.1 )
            {
                fp("%g %g\n", P, param.sig_DelP(P));
            }
        }
        
        
        return 0;
    }
    catch(...)
    {
        std::cerr << "Error !" << std::endl;
    }
    
    return -1;
}